On : “ Numerical integration of related Hankel trans - forms by quadrature and continued fraction expansion ”

Chave presents an excellent algorithm and computer subprogram that evaluate Hankel transforms by quadrature and continued fraction summation. This subprogram would be a valuable addition to any mathematical computer library. Chave refers to digital filter (convolution) methods as “the standard numerical approach to the computation of Hankel transforms,” and makes numerous references to my published work (Anderson, 1979, 1982). However, some readers may misinterpret some of his statements regarding convolution methods, which I hope to clarify by this discussion. On p. 1671, Chave states (referring to Anderson, 1979) that “Reasonable (5-figure) accuracy is typically achieved for monotonic, rapidly decreasing kernel functions at moderate values of p.” This statement implies a stronger condition imposed on the kernels than defined in Anderson (1979, p. 1289), where only continuous, bounded kernels are required. Rapid convergence of the convolution sum is achieved if the kernel function is also decreasing (Anderson, 1982) but the rate of decrease is usually immaterial, inasmuch as extremely rapid decaying digital filter responses are used. The 5-figure accuracy is typical of singleprecision (32 or 36-bit floating-point) implementations as exemplified in subprograms ZHANKS (Anderson, 1979) and HANKEL (Anderson, 1982), where the best relative errors are approximately 10m6. Chave’s brief reference on p. 1671 to “adaptive and lagged convolution” in Anderson (1982) failed to mention that a double-precision (64-bit floating-point) related and lagged con~d~tion subprogram (DHANKL) is available, where best relative errors are approximately IO“. He further states on p. 1671: “For some types of problems the digital filter method is less useful; examples occur at very small values of the range p, ...> and when high numerical precision is required.” With either HANKEL or DHANKEL, very small values of p 2 0 can be accommodated by lagged convolution as discussed by Anderson (1982, p. 366: “Proceeding to the limit”). In some applications the need for double-precision and very small arguments may be avoided by using suitable transformations of the Hankel integral. For example, the electromagnetic (EM) induction problem discussed by Anderson (1979, p. 12921293) used a normalized induction number (B = p/6, where 6 is the skin depth of the medium) instead of the usual distance p parameter in the Bessel argument of the Hankel transform. This approach usually avoids small transform B arguments, where B is typically in the range [.Ol, lo] for the quasi-static assumption. Algebraically divergent Hankel transforms of the type discussed by Anderson (1979, p. 1293) and Frischknecht (1967, p. 67) can be replaced with rapidly convergent ones by subtracting a known homogeneous half-space term under the integral and adding an equivalent analytic expression outside the integral. For many practical EM problems, combining these two substitutions will result in Hankel transforms that converge rapidly for a moderate induction number range; they are therefore computationally tractable using the singleprecision convolution subprograms ZHANKS, HANKEL, or similar routines. The oscillatory bounded kernels mentioned on p. 1674 can be transformed successfully with moderate to large arguments in ZHANKS, HANKEL, or DHANKL, if the highest frequency content of the kernel function does not approach or exceed the digital filter’s sampling or Nyquist frequency (see Anderson, 1982, p. 346, p. 362 regarding oscillating functions). Highly oscillatory kernels are rare in practice, but when they do occur, the quadrature algorithm Chave proposed would generally be more useful and accurate than convolution methods. On p. 1674, the statement, “Only the first pair of integrals (6H7) can be handled by the digital filter method” is misleading in view of the above discussion. Integrals (8)(9) and (12t (13) can also be numerically transformed using convolution methods, since the convolution input (kernel) functions meet the basic requirements defined in Anderson (1982, p. 346) and Anderson (1979, p. 1289). Equations (lO)(ll) have increasing (unbounded) kernels, and would certainly fail for any argument value if used directly in ZHANKS or DHANKL. Algebraically divergent integrals can often be converted to convergent form using transformations as previously mentioned; however, I will not consider (1 OH 11) any further. To illustrate the above discussion, subprogram DHANKL (Anderson, 1982) was run on a VAX-l l/780 VMS3.5 system in double-precision complex arithmetic (COMPLEX*16) for the same examples presented in Chave’s Table 2. The divergent types (lO)( 11) were excluded, as well as the oscillating types (12)(13) for argument R = 0.05. My DHANKL results are listed in Table I. The tolerance factor (TOL) used in Table I was selected the same as parameter RERR in Chave’s Table 2. The number of kernel function evaluations used in DHANKL for direct convolution is denoted by NF in Table I, and NF = 0 denotes related convolution was used. Note that the accuracy in NUMERR for oscillatory integrals NN = 7 and NN = 8 at R = 2.0 is only good to about three figures; however, these same integrals at R = 100.0 are nearly equivalent to those in Chave’s Table 2. All remaining results in Table I clearly show comparable accuracy with Chave’s Table 2 and with respect to the computed exact values and requested accuracy (see TOL in Anderson, 1982, p. 352 and double-precision version, p. 364). It is stated on p. 1674 that “Table 2 shows similar computations [as in Table I] with RERR = lolo, a far more stringent requirement with a concomitant increase in compu-